############
# packages #
############
library(PanelMatch)
library(plm)
library(dplyr)
library(ggplot2)
library(patchwork)
library(labelled)

###############
# set-up data #
###############
# read dataset
dat.anlsys <- read.csv("Party data for analysis.csv")
dat.anlsys <- as.data.frame(dat.anlsys)
# turn variables from integer to numeric
dat.anlsys <- dat.anlsys %>% mutate_if(is.integer, as.numeric)
# keep time and unit id as integer
dat.anlsys$period <- as.integer(dat.anlsys$period)
dat.anlsys$v2paid <- as.integer(dat.anlsys$v2paid)

######################
# Matching procedure #
######################
### perform matching with refinement
## moderate inclusiveness measure
# generate panel data for PanelMatch
dat.anlsys.dlgts <- PanelData(panel.data = dat.anlsys, time.id = "period",
                              unit.id = "v2paid",outcome = "v2pavote",
                              treatment = "prmrs.dlgts")
# perform matching 
PM.results.dlgts <- PanelMatch(panel.data = dat.anlsys.dlgts, lag = 2, match.missing = TRUE, 
                               refinement.method = "ps.weight",
                               covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                 v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                 I(lag(gov.mmbr,1)) + I(lag(prmrs.dlgts.swtch.other,1)) + 
                                 e_gdppc + v2paallian, 
                               qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)

## strict inclusiveness measure
# generate panel data for PanelMatch
dat.anlsys.mmbrs <- PanelData(panel.data = dat.anlsys, time.id = "period",
                              unit.id = "v2paid",outcome = "v2pavote",
                              treatment = "prmrs.mmbrs")
# perform matching 
PM.results.mmbrs <- PanelMatch(panel.data = dat.anlsys.mmbrs, lag = 2, match.missing = TRUE, 
                               refinement.method = "ps.weight",
                               covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                 v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                 I(lag(gov.mmbr,1)) + I(lag(prmrs.mmbrs.swtch.other,1)) + 
                                 e_gdppc + v2paallian, 
                               qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)

#########################
# estimation (Figure 2) #
#########################
## set seed for boostrap se reproducability
set.seed(1234)
## calculate estimates (must run matching code first for PanelMatch objects)
PE.results.dlgts.90 <- PanelEstimate(sets = PM.results.dlgts, panel.data = dat.anlsys.dlgts, confidence.level = 0.90)
PE.results.dlgts.95 <- PanelEstimate(sets = PM.results.dlgts, panel.data = dat.anlsys.dlgts, confidence.level = 0.95)
PE.results.mmbrs.90 <- PanelEstimate(sets = PM.results.mmbrs, panel.data = dat.anlsys.mmbrs, confidence.level = 0.90)
PE.results.mmbrs.95 <- PanelEstimate(sets = PM.results.mmbrs, panel.data = dat.anlsys.mmbrs, confidence.level = 0.95)
## present numeric outcomes
summary(PE.results.dlgts.90, verbose = FALSE)
summary(PE.results.dlgts.95, verbose = FALSE)
summary(PE.results.mmbrs.90, verbose = FALSE)
summary(PE.results.mmbrs.95, verbose = FALSE)
## coefficient plots for both measures
# extract data frame from estimate
out.dlgts.95 <- summary(PE.results.dlgts.95, verbose = FALSE)
out.dlgts.90 <- summary(PE.results.dlgts.90, verbose = FALSE)
out.mmbrs.95 <- summary(PE.results.mmbrs.95, verbose = FALSE)
out.mmbrs.90 <- summary(PE.results.mmbrs.90, verbose = FALSE)
model <- c("dlgts","dlgts","mmbrs","mmbrs")
att <- c(out.dlgts.95[1], out.dlgts.90[1], out.mmbrs.95[1], out.mmbrs.95[1])
ci.upr <- c(out.dlgts.95[3], out.dlgts.90[3], out.mmbrs.95[3], out.mmbrs.90[3])
ci.lwr <- c(out.dlgts.95[4], out.dlgts.90[4], out.mmbrs.95[4], out.mmbrs.90[4])
ci.lvl <- c(95,90,95,90)
dat.PEresults.plot <- data.frame(att,model,ci.upr,ci.lwr,ci.lvl)
## plot results
PEresults.plot <- ggplot(dat.PEresults.plot, aes(y=att,x=model,group=model)) +
  geom_point(size=2.5) +
  geom_linerange(dat=filter(dat.PEresults.plot,ci.lvl==95),aes(ymin=ci.lwr,ymax=ci.upr),linewidth=0.5) +
  geom_linerange(dat=filter(dat.PEresults.plot,ci.lvl==90),aes(ymin=ci.lwr,ymax=ci.upr),linewidth=1) +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_y_continuous(limits=c(-13,13), breaks=seq(-10,10,5)) +
  scale_x_discrete(labels=c("Moderate", "Strict")) +
  labs(y="Treatment Effect on Vote Share", x="Inclusiveness Measure") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title = element_text(size=12),
    axis.text = element_text(size=11),
  )
## save plot (Figure 2)
ggsave(file="Figure2.pdf", plot=PEresults.plot, width=5, height=3.5)
